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Studies of random close packing of spheres have advanced our knowledge about the structure 
of systems such as liquids, glasses, emulsions, granular media, and amorphous solids. When these 
systems are confined their structural properties change. To understand these changes we study 
random close packing in finite-sized confined systems, in both two and three dimensions. Each 
packing consists of a 50-50 binary mixture with particle size ratio 1.4. The presence of confining walls 
significantly lowers the overall maximum area fraction (or volume fraction in three dimensions). A 
simple model is presented which quantifies the reduction in packing due to wall-induced structure. 
This wall-induced structure decays rapidly away from the wall, with characteristic length scales 
comparable to the small particle diameter. 



I. INTRODUCTION 

In 1611 Kepler conjectured that the most efficient 
packing of spheres is the face center cubic packing {fee), 
with packing fraction (f>f cc = tt/^/18 pQ. In 1831, Gauss 
provided a partial proof to this conjecture, and more re- 
cently Hales presented a more complete proof [5] that 
is still being validated pQ. However, the fee packing 
is highly ordered, and in many cases random packings 
are of interest, as they are easy to create. For exam- 
ple, a large range of packing fractions have been found 
for granular particles with a minimum mechanically sta- 
ble volume fraction (j> r i p w 0.55 [3j|4j|5], termed random 
loose packing. The maximum volume fraction for a ran- 
domly packed 3D system is 4> rcp ~ 0.64, termed random 
close packing {rep) [3j IH El U\ [8] . 

The earliest known experimental study on the density 
of rep was performed by Bernal and Mason [5] . In their 
experiment they repeatedly shook and compressed a rub- 
ber balloon full of spheres for sufficiently long enough 
time to reach a very dense state at (f> rcp — 0.637. This 
result is sensitive to the preparation method, although 
the final volume fractions found are always close to 0.64 
[9]. Computational studies are especially sensitive to pro- 
tocol, with (j> rcp ranging between 0.64 and 0.68 [10| lllj. 
A broader range of mechanically stable packing fractions 
can be obtained by consider packings consisting of a 
polydisperse mixture of spheres or packings with non- 
spherical particles [12l HH HH [HI [T6l [17l UHl [HI [2Q]. 
Some theoretical progress explaining random close pack- 
ing has occurred [HJ H21 [21], although this is far 
from complete [TJ. All of these studies are relevant to 
a wide range of problems including the structure of liv- 
ing cells [5], liquids [BJ, granular media [STJ [52], emul- 
sions [22, glasses [25], and amorphous solids [27] , 

While all of the above studies focus on infinite systems, 
real systems have boundaries and often these bound- 
aries are important. Furthermore, in many cases sam- 
ples are confined between closely spaced boundaries, with 
the confining size being only a few characteristic parti- 
cle sizes across 28J. For example, when a liquid is con- 
fined its structure is dramatically changed with particles 



forming layers near the wall, which ultimately affects the 
properties of the liquid [29] Hfl EJ [32] [33l [34]- The 
shearing of confined dense colloidal suspensions shows 
the emergence of new structures not seen before [35] . The 
flow of granular media through hoppers [36] [37] or sus- 
pensions through constricted micro- and nanofluidic de- 
vices [351 [Ml HOI SI] can jam and clog, costing time and 
money. Some studies examined the packing of granu- 
lar particles in narrow silos, focusing on stresses between 
particles and the walls [HJ [331 S31 H5J- Other studies 
noted that the particle packing in silos is layered near 
the walls [151 H7] . However, the influence of boundaries 
on (j) rcp has not been studied previously. 

We present computer simulated rep packings in con- 
fined geometries. In particular, we study binary mix- 
tures to prevent wall-induced crystallization |48, 4"9"1I50] . 
We create two-dimensional (2D) and three-dimensional 
(3D) packings between two parallel walls, with periodic 
boundary conditions in the other directions. Confine- 
ment significantly modifies the rep states, with lowered 



values for 



J rcp 



reflecting an inefficient packing near the 



walls. This inefficient packing persists several particle 
diameters away from the wall, although its dominant ef- 
fects are only within 1-2 diameters. 

These results will be useful for understanding many 
other confined systems. For example, many experiments 
study how confinement modifies the glass transition; 
samples which have a well-characterized glass transition 
in large samples show markedly different properties when 
confined to small samples [251 |5Dl EH [53 [53J [5U [531 [Sg 
157] , However, it's not clear if these changes are due to 
boundaries or finite size effects [58 . Our results show 
that boundaries significantly modify the packing, which 
may in turn modify behavior of these confined molecular 
systems [2"5] . 

The manuscript is organized as follows. Section [TT] out- 
lines the algorithm we use to generate confined rep states. 
Section |III| shows how the total packing fraction, parti- 
cle number density, and local order of confined rep states 
change with confining thickness and distance from the 
confining boundary. Finally, Sec. IV provides a simple 
model that predicts the packing fraction dependence with 
confinement. 
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FIG. 1: A flow chart outlining our algorithm for computing 
rep configurations. 



II. METHOD 

Our aim is to quantify how a confining boundary alters 
the structure of randomly closed packed (rep) disks in 
2D and spheres in 3D, and in particular to study how 
this depends on the narrowest dimension. This section 
presents our algorithm for 2D packings first, and then 
briefly discusses differences for the 3D algorithm. 

In 2D, our system consists of a binary mixture of disks 
containing an equal number N/2 large disks of diameter 
di and small disks of diameter d s with size ratio a — 
di/d s — 1.4. For each configuration, disks are packed into 
a box of dimensions L x by L y , with a periodic boundary 
along the x-direction and a fixed hard boundary (a wall) 
along the y-direction. 

Each configuration is generated using a method 
adapted from Xu et al. [8] which is an extension of a 
method proposed by Clarke and Wiley jSH] ■ This method 
is briefly summarized in Fig.[T] Infinitesimal particles are 
placed in the system, gradually expanded and moved at 
each step to prevent particles from overlapping. When a 
final state is found such that particles can no longer be 
expanded without necessitating overlap, the algorithm 
terminates. Near the conclusion of the algorithm, we 
alternate between expansion and contraction steps to ac- 
curately determine the rep state. 

In particular, while the final state found is consistent 
with hard particles (no overlaps allowed), the algorithm 
uses a soft potential at intermediate steps [8] , given by 



V ( r ij) = o I 1 ~ r ij/dij) 6 (1 - Tij/dij) , 



(1) 



disk i and j, e is a characteristic energy scale (e = 1 for 
our simulations), dij — (di + dj)/2, and 6 (1 — rij/dij) is 
the Heaviside function making V nonzero for r^j < dij. 
Simulations begin by randomly placing disks within a box 
of desired dimensions and boundary conditions with the 
initial diameters chosen such that 4>initiai 4>rcp- In the 
initial state particles do not overlap and the total energy 
E = 0. 

Next all disk diameters are slowly expanded subject to 
the fixed size ratio a = 1.4 and <p changing by 8<p per 
iteration; we start with 6<j) = 10~ 3 . After each expansion 
step, we check if any disks overlap, by checking the condi- 
tion 1 — rij/dij > e r — 1CP 5 for each particle pair. Below 
this limit, we assume the overlap is negligible. If any par- 
ticles do overlap (E > 0), we use the non- linear conjugate 
gradient method [60] to decrease the total energy by ad- 
justing the position of disks so they no longer overlap 
(E = 0). In practice, one energy minimization step does 
not guarantee we have reached a minimum within the de- 
sired numerical precision. Thus this step can be repeated 
to further reduce the energy if E > 0. We judge that we 
have reached a nonzero local minimum if the condition 
\\VE\\/(2N) < e E = 1CT 7 is found, where ||V£|| is the 
magnitude of the gradient of E. Physically speaking this 
is the average force per particle, and the threshold value 
(10~ 7 ) leads to consistent results. 

If we have such a state with E > 0, this is not an 
rep state as particles overlap. Thus we switch and now 
slowly contract the particles until we find a state where 
particles again no longer overlap (within the allowed tol- 
erance). At that point, we once again begin expansion. 
Each time we switch between expansion and contraction, 
we decrease 5<f) by a factor of 2. Thus, these alternating 
cycles allow us to find an rep state of non-overlapping 
particles (within the specified tolerance) and determine 
4> rC p to high accuracy. We terminate our algorithm when 
5(f) < <5</> m i n = 1CP 6 . In practice, we have tested a variety 
of values for the thresholds e r ,es, and <50 m i n and find 
that our values guarantee reproducible results as well as 
reasonably fast computations. Our algorithm gave an 
average packing fraction of 4> rcp — 0.8420 ± 0.0005 for 
40 simulated rep states containing 10,000 particles with 
periodic boundary conditions along both directions. Our 
value of (p r cp is in agreement with that found by Xu et 
al. 0. 

The above procedure is essentially the same as Ref . [H] ; 
we modify this to include the influence of the boundaries. 
To add in the wall, we create image particles reflected 
about the position of the wall; thus particles interact 
with the wall using the same potential, Eqn. [T] 

Additionally, we wish to generate packings with pre- 
specified values for the final confining height h = L y /d s . 
(This allows us to create multiple rep configurations with 
the same h.) We impose h by affinely scaling the system 
after each step, so that the upper boundary is adjusted 
by L y — hd s and each disk's y-coordinate is multiplied by 



the ratio L 



/L 



y,l1 



where L yi and L y i+ i 



are the con- 



where nj is the center to center distance between two fining widths between two consecutive iterations. Thus 
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FIG. 2: (Color online). Illustrations of 2D and 3D config- 
urations generated using the algorithm described in Sec. [TT] 
(a) 2D configuration for h = 10. (b) 2D configuration with 
h — 20. (c) 3D configuration with h — 5 where blue (dark 
gray) represents big particles and green (light gray) represents 
small particles, (d) 3D configuration with h = 10. 



while d s gradually increases over the course of the simu- 
lation, L y increases proportionally so that the nondimen- 
sional ratio h is specified and constant. Some examples 
of our final rep states are shown in Fig. [2j 

To ensure we will have no finite size effects in the peri- 
odic direction, we examined <j) rC p for different h and L x , 
and found <fi rcp (h) to be independent of L x for 3 < h < 30 
if L x /d s > 40. Thus we have chosen N for each simula- 
tion to guarantee L x /d s ~ 50. 

In 3D, our system consists of a binary mixture of 
spheres containing an equal number N/2 large spheres 
of diameter di and small spheres of diameter d s with a 
size ratio a = di/d s = 1.4. Spheres are packed into a box 
of dimensions L x by L y by L z , with periodic boundaries 
along the x- and z-directions and a fixed hard boundary 
along the y-direction. Each configuration is generated us- 
ing the same particle expansion and contraction method 
described above and the same initial values for 84> and the 
terminating conditions. For each configuration L x = L z , 
h = Ly/d s , and N is chosen so that L x /d s > 10. Our 
choice of L x /d s > 10 is not large enough to avoid finite 
effects. However, in order to acquire the large amount 
of data needed in a reasonable amount of time we in- 
tentionally choose a value of L x /d s below the finite size 
threshold. Trends observed in the 2D analysis will be 
used to support that any similar trends seen in 3D are 
real and not the result of the finite periodic dimensions. 
Note that in 3D we will show cases where h > L x /d s 
resulting in the confining direction being larger than the 
periodic direction, and this may affect the structure of fi- 
nal configurations; however, we will not draw significant 
conclusions from those data. 

Overall, it is not known if this algorithm produces 
mathematically rigorously defined random close packed 
states [5J H3 |BTJ [55]. However, the goal of this article is 
to determine empirically the properties of close-packed 
states in confinement, and we are not attempting to ex- 
tract mathematically rigorous results. For example, we 
are not as interested in the specific numerical values of 
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FIG. 3: (Color online). The black curve is the average packing 
fraction (f> found by averaging at least 10 2D configurations 
together for various confining widths h; recall that h has been 
nondimensionalized by d s , the small particle diameter. The 
red curve (dark gray) is a fit using Eqn.[5]which finds (/> rcp = 
0.842 in the limit h — > oo; the value for 4> rcp is indicated 
by the black dashed line. The green (light gray) data points 
are <j)(h) computed for many configurations with the confining 
wall replaced by a periodic boundary. The inset is a magnified 
view of the region for ft < 6 to better show the large variations 
within this range. The vertical lines in the inset are located 
at "special" h values where peaks and plateaus appear. 



4>rcp that we obtain, but rather the qualitative depen- 
dence on h. As noted in the introduction, different com- 
putational and experimental methods for creating rep 
systems have different outcomes, and so it is our qual- 
itative results we expect will have the most relevance. 

Note that for the remainder of this paper, we will drop 
the subscript rep, and it should be understood that dis- 
cussions of (f> refer to the final state found in each simu- 
lation run, <p rcp (h). 



III. RESULTS 

A. 2D Systems 

We begin by generating many 2D configurations with 
h between 3-30 and computing the packing fraction for 
each, as shown by the black curve in Fig. [3] This plot 
shows that confinement lowers <p, with the influence of 
the walls being increasingly important at lower h. The 
lowering of <p with confinement is most likely due struc- 
tural changes in the packing near the confining boundary. 
We know that any alteration in particle structure from 
a rep state must be "near" the wall because as h — > oo 
we expect to recover a packing fraction of </> rcp , implying 
that in the infinite system the "middle" of the sample 
is composed of an rep region. Extrapolating the data in 
Fig. [3] to h — > oo we find 4>h^oo = <t>rcp = 0.842 which is 
essentially a test of our method. The extrapolation (red 
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FIG. 4: A plot of the number density p(y) for 100 2D config- 
urations at h — 30 averaged together. The plot is constructed 
by treating the small and big particles separately and using 
bins along the confining direction of width 8y = 0.1d s . 



curve in Fig. [3| was carried out by assuming that to first 
order </> ~ 4>h^oo — C /h for large h, where </>/i_>oo = 4>rcp 
(the bulk value for the rep packing) and C is a fitting 
parameter. 

The data in Fig.wmegin to deviate from the fit for h < 
6, and furthermore (p(/i) is not monotonic. While some of 
the variability is simply noise due to the finite number of 
disks N used in each simulation, some of the variability 
is real. The inset in Fig. [3] shows a magnified view of 
the region 3 < h < 6. The vertical lines in this inset are 
located at specific values of h that can be expressed as the 
integer sums of the two particle diameters. For instance, 
the first vertical line near the y-axis is located at h = 
2d a + di . These lines are placed at some h values where 
(f>(h) has notable spikes or plateaus. These lines suggest 
that there exist special values of h where the confining 
thickness is the right width so that particles can pack 
either much more efficiently or much less efficiently than 
nearby values of h. fntriguingly, these special h values 
are not the set of all possible integer sums, but instead 
only the selected few drawn in the figure. For example, 
there is no apparent feature near h = 4 (recall that h is 
nondimensionalized by d s ). Somewhat surprisingly, the 
peaks correspond to integer combinations of d s and di, 
rather than combinations such as \/3/2d s . The latter 
would suggest hexagonal packing, the easiest packing of 
monodisperse disks in 2D; whereas the observed peaks of 
4>(h) suggest square-like packing. 

To measure structural changes in particle packing as 
a result of confinement we start by examining the varia- 
tions in the local number density p with distance y from 
the confining wall. Figure [4] is a plot of p(y) for 100 con- 
figurations averaged together at h — 30. This plot shows 
oscillations in particle density which decay to a plateau. 
The oscillations near the wall are indicative of particles 
layering in bands. Above y > 6d s , noise masks these 



oscillations. This supports our interpretation, that con- 
finement modifies the structure near the walls but not 
in the interior. Furthermore, the rapidity of the decay 
to the plateau seen in Fig. [4] suggests that confinement 
is only a slight perturbation to systems with overall size 
h>6. 

The details of the density profiles in Fig. [4] also sug- 
gest how particles pack near the wall. The small particle 
density (solid line) has an initial peak at y = 0.5d s , indi- 
cating a large amount of small particles in contact with 
the wall, as their centers are one radius away from y = Q. 
Likewise, the large particle density (dashed line) has its 
initial peak at y = 0.7d s = 0.5d/, indicating that those 
particles are also in contact with the wall. This is con- 
sistent with the pictures shown in Fig. [2][a,b), where it is 
clear that particles pack closely against the walls. Exam- 
ining again the small particle number density in Fig. [4] 
(solid line), the secondary peaks occur at y — l.5d s and 
y = l.9d s = 0.5d s + l.0di, which is to say either one small 
particle diameter or one large particle diameter further 
away from the first density peak at y = 0.5d s . This again 
is consistent with particles packing diamctcr-to-diamctcr, 
rather than "nesting" into hexagonally packed regions. 
Similar results are seen for the large particles (dashed 
line) which have secondary peaks at y = l.0d s + 0.5d; 
and y — 2.ld s — l.bdi. 

To confirm that these density profile results apply for 
a variety of thicknesses h, and more importantly to see 
how these results are modified for very small h, we use 
an image representation shown in Fig. [5| To create this 
image, density distributions of different h are each sep- 
arately rescaled to a maximum value of 1. Every data 
point within each distribution is then made into a gray 
scale pixel indicating its relative value; black is a relative 
value of 1, and white is a relative value of 0. The vertical 
axis is the confining width and the horizontal axis is the 
distance y from the bottom wall. Each horizontal slice 
(constant h) is essentially the same sort of distribution 
shown in Fig. [4] The white space on the right side of the 
figure arises because the distribution is only plotted for 
the range < y < h/2. The distributions are symmetric 
about y = h/2, and by averaging the distribution found 
for the range < y < h/2 with the distribution found for 
the range h/2 < y < h, the statistics are doubled. The 
area shown in the box is a magnified view of that region 
where the full range < y < h is being shown. 

In Fig. [5] there are vertical strips of dark areas, once 
again indicating that particles arc forming layers. The 
width of these strips widens and the intensity lessens far- 
ther from the wall. In each plot, the first vertical black 
strip is sharply defined and located at one particle radius, 
illustrating that small and big particles are in contact 
with the wall. Finally, the location and width of each 
layer remains essentially the same for different h, sug- 
gesting that layering arises from a constraint imposed by 
the closest boundary. Given that the first layer of parti- 
cles always packs against the wall, this imposes a further 
constraint on how particles pack in the nearby vicinity. 
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FIG. 5: An image representation constructed for the purpose 
of comparing 2D p(y) distributions at many different h. The 
intensities have been logarithmically scaled. The vertical pixel 
width is 0.1 and for the left plot the horizontal pixel width is 
0.2 and for the right plot the horizontal pixel width is 0.14. 



The consistency in the location and width of the sec- 
ond layer for all h demonstrates that the constraint of 
the first layer always produces a similar packing in the 
second layer, essentially independent of h. Continuing 
this argument, each layer imposes a weaker constraint on 
the formation of a successive layer, allowing for the local 
packing to approach rep far from the wall. 

In the magnified views of Fig. [5j the vertical dark lines 
show the layering of particles induced by the left bound- 
ary and the angled dark lines show the layering of par- 
ticles induced by the right boundary. We see that for 
small h these sets of lines overlap and intersect, meaning 
that there is a strong influence from one boundary on the 
packing within the layers produced by the other bound- 
ary. This may explain the variations seen in <j>(h) for 
small h in Fig. [3] In particular, it is clear that at certain 
values of h, the layers due to one wall are coincident with 
the layers due to the other wall, and this suggests why 
4>(h) has a higher value for that particular h. Given that 
the layer spacings correspond to integer combinations of 
d a and di, the coincidence of layers from both walls will 
correspond to integer combinations of d s and di , and this 
thus gives insight into the peak positions shown in the 
inset of Fig. [3j 

As described above, the influence of the walls dimin- 
ishes rapidly with distance y away from the wall. In par- 
ticular, for the local number density p{y), we observe that 
the asymptotic limit p(y — > oo) = 0.362 for the curves 
shown in Fig.[4]is in agreement with the theoretical num- 
ber density of an rep configuration p rcp — 4</> rcp /7r(l + cr). 
To quantify the approach to the asymptotic limit, we de- 
fine a length scale from a spatially varying function f(y) 
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FIG. 6: (Color online). Drawings illustrating the conceptual 
meaning of (a) (b) i/jq, and (c) t[>7. Darker colored particles 
have neighbors that are packed more like an ideal regular n- 
sided polygon as compared to lighter drawn particles. The 
configuration of particles is the same for all panels, and are 
drawn from a simulation with h = 10. Note in (b) that there 
are no large patches of high tpe , demonstrating that there are 
no large crystalline domains. 



using: 



A = 



fy[f(y) - f(y -> °°)] 2 dy 

J[f(y)-f(y^^)] 2 dy ' 



(2) 



In this equation f(y) is an arbitrary function where the 
value of A quantifies the weighting of f(y). For sim- 
ple exponential decay f(y) — Ae^ v ^ x \ Eqn. [2] gives 
A = A'/2. Using f(y) = p(y) we find A = 0.854 and 
A = 0.72d s for the small particle curve and big particle 
curve in Fig. [4] respectively, suggesting that the transi- 
tion from wall-influenced behavior to bulk rep packing 
happens extremely rapidly. 

To further investigate the convergence of the local 
packing to rep more closely we analyze the local bond or- 
der parameters ?p n , which for a disk with center of mass 
Ts are defined as 



J_\p e mS(r l3 ) 

3 



(3) 



The sum is taken over all j particles that are neighbors 
of the ith particle, 0(ry) is the angle between the bond 
connecting particles i and j and an arbitrary fixed refer- 
ence axis, and rib is the total number of i - j bonds |63j . 
The magnitude of ip n is bounded between zero and one; 
the closer the magnitude of ip n is to 1, the closer the lo- 
cal arrangement of neighboring particles are to an ideal 
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FIG. 7: (Color online), (a) is a plot of the local num- 
ber density p(y) for 2D configurations of big and small par- 
ticles separately, (b) - (d) are plots of {^n){y) for small 
(green/light gray) and big particles (blue/dark gray) sepa- 
rately, and both sizes together (light purple/medium gray) 
where (b) is {ips}, (c) is (tpe), and (d) is (V>7). The length 
scales determined from these curves for small, large, and both 
species are A 5 . s = 1.2,A 5ii = 1.1, A 5 . b = 1.4, A6, s = 0.8,A 6 ,; = 
0.9, A 6 , 6 = 0.8,A 7 , S = 1.1, A 7 j = 0.7, and A 7 ,i, = 1.0 (all in 
terms of d s ). 



n-sided polygon. Figures [6|a-c) are drawings illustrating 
the concept of ip n using a 2D configuration with h = 10. 
Particles with larger tp n are drawn darker. These figures 
have no large clusters of dark colored particles, demon- 
strating that there are no large crystalline domains (i.e. 
particles are randomly packed). 

For a highly ordered monodisperse packing (ip e ) would 
be the most appropriate choice for measuring order be- 
cause of the ability for monodisperse packings to form 
hexagonal packing. However for a bidisperse packing 
with size ratio a — 1.4, the average number of neighbors 
a small particle will have is 5.5 and the average number 
of neighbors big particles will have is 6.5. Therefore, a 
bidisperse packing of this kind will have a propensity to 
form local pentagonal, hexagonal, and heptagonal pack- 
ing, and to properly investigate how the local packing 
varies we examine (ips), (ip6), and ("07). We compute 
the average values (i/'s), (ipe), and (1P7) for all configura- 
tions as a function of y, and averaging together all (ipn) 
distributions for configurations with h > 16 to improve 
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FIG. 8: The black data points are the average packing frac- 
tions of 3D configurations at various h. The red (dark gray) 
curve is a fit from our model Eqn. [5] For each h at least 10 
configurations were averaged together. 



statistics. This averaging can be justified by considering 
that oscillations in p(y) in Fig. [4] for y/d s > 10 are quite 
small. Thus this averaging improves our statistics for the 
range < y/d s < 5 where the largest oscillations occur, 
without skewing the data. In the end nearly 10,000 con- 
figurations are averaged together, producing the curves 
shown in Fig.mb-d). This figure shows the spatial vari- 
ations of (ip5}7\tpe}i and (^7) for small and big particles 
separately and both particles combined. All curves show 
fluctuations that decay with distance from the wall, and 
show local order within and between layers. Figure [7^ a) 
has been added so comparisons between the location of 
the oscillations in p(y) and (ip n )(y) can be made. 

For the most part, each successive layer has less orien- 
tational order than the previous layer with (ip n ) eventu- 
ally decaying to an asymptotic limit. To characterize a 
length scale for these curves we compute A using Eqn. [2] 
for each curve shown in Fig[7]jb-d). From the nine curves, 
we find that the mean value of A — (1.00 ± 0.24)d s . The 
length scales found for these curves are once again less 
than the largest particle diameter. No striking differ- 
ence is found between the different order parameters or 
between the different particle sizes; specific values of A 
are given in the figure caption. (Note that the asymp- 
totic limits of all (ip n ) plots are in agreement with the 
average values found for 40 unconfined 10,000 particle 
simulations averaged together, confirming that the local 
packing converges to an rep arrangement far from the 
walls.) 

Next, we wish to distinguish the structural influence of 
the flat wall from the finite size effects. We perform simu- 
lations where the confining wall is replaced by a periodic 
boundary with periodicity h; thus particles cannot form 
layers. In this case, the packing fraction still decreases 
as h is decreased, as shown by the green curve (light 
gray) in Fig. [3j although the effect is less striking than 
for the case with walls (black curve) . A likely explanation 
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for the decrease in <f> with confinement is the long range 
structural correlations imposed along the constricted di- 
rection; in other words, if there is a particle located at 
(x, y) that particle is mirrored at (x,y — h) and (x,y + h) 
by the periodic boundary. We know from the pair corre- 
lation function [12] E] of rep configurations that struc- 
tural correlations exist over distances of many particle 
diameters, although of course these are weak at larger 
distances. Thus the periodicity forces a deviation from 
the ideal rep packing, that becomes more significant as h 
decreases. By definition rep is the most random densely 
packed state, and thus any perturbations away from this 
state must have a lower packing fractions. However, this 
is not nearly as significant as the constraint imposed by 
the flat wall, as is clear comparing the green (light gray) 
data and the black data in Fig. [3] 



B. 3D Systems 

We next consider 3D confined systems. We start by 
investigating <j)(h), shown as the black points in Fig. [8] 
As observed in the 2D case, 4> is reduced as a result of 
confinement. However, unlike the 2D system, there does 
not appear to be a scries of "special values" of h that give 
rise to peaks and plateaus, other than a hump near h — 
3.75. The lack of substructure may be due to the smaller 
size in x and z, in contrast with the 2D simulations which 
had large sizes in the unconfined direction. 

Next we investigate the local number density p(y) for 
h = 25, shown in Fig. |9ja). The data are constructed by 
averaging together 100 configurations. The curve shows 
fluctuations that decay with distance from the wall, even- 
tually reaching a plateau. Using Eqn. [2] we obtain de- 
cay lengths A3B = 0.77d s and 0.73<i s for the small and 
large particle curves respectively. These length scales 
are similar to the length scales obtained in the 2D case 
(A2D = 0.85<i s and 0.72d s for small and large particles). 

To compare all 3D p(y) distributions for different h we 
construct the image representation used to compare 2D 
configurations in Fig. [5} T he data for the 3D configura- 
tions are shown in Fig. |10| Again there are dark vertical 
strips arising from particles forming layers near the wall. 
Like in 2D, the density approaches the "bulk" rep value 
far from the wall. 

In 2D, we also noted that the structure is modified 
near the wall, as measured by the ipn order parameters. 
To investigate structural ordering in 3D, we use a local 
structural parameter sensitive to ordering [53] [55] . We 
start by defining 



?i,6 



— y 

J 3 



(4) 



In the above equation m — {—6, 0, 6}, and thus qi$ 
is a 13 element complex vector which is assigned to every 
particle i in the system. The sum in Eqn. [4] is taken over 
the j nearest neighbors of the ith particle, rij is the total 




i> 1.4 - 



yld s 



FIG. 9: (a) is plot of p(y) for 3D configurations for small and 
big particles separately. The plot is constructed using bins of 
width 8y = 0.1d s along the confining directions, (b) is a plot 
of the average number of ordered bonds (Nb)(y). 



number of neighbors, and K is a, normalization constant 
so that qifi ■ cji.e — 1. For two particles i and j that 
are nearest neighbors, YQ m (0ij , <fiij) is the spherical har- 
monic associated with the vector pointing from particle 
i to particle j, using the angles 0y and 4>ij of this vector 
relative to a fixed axis. Next, any two particles m and 
n are considered "ordered neighbors" if q m ,6 ' Qn, 6 > 0.5 
[641 165) . Finally, we quantify the local order within the 
system by the number of ordered neighbors Nb a particle 
has. 



Figure |9|b) is a plot of the average number of ordered 
neighbors particles have (Nf,) as a function of distance 
y from the wall. In comparison with Fig. 9ja), this plot 
shows that local order is mostly seen within layers, not 
between layers. Also we see that (Nb) converges to an 
asymptotic value of ~ 1.3, confirming that the system 
is disordered (values of Nb > 8 are considered crystalline 
[65]). We use Eqn. [2] to characterize a length scale for the 
decay in (Nb), giving A = l.3d s . The asymptotic limit of 
(Nb)(y) in Fig. ^b) agrees with the average value of Nb 
found for 15 large simulations with 2,500 particles and 
periodic boundary conditions, confirming that the local 
structure in the confined case converges to the bulk rep 
state far from the walls. 

Our results show that in both 2D and 3D, confinement 
induces changes in structural quantities near the walls, 
with a decay towards the "bulk" values characterized by 
length scales no larger than di. The only prior work 
we are aware of with related results are a computational 
study [46] and an experimental study [47] of collections 
of monodisperse particles confined in a large silo. The 
simulation by Landry et al. primarily focused on the force 
network within the silo. They show one plot of the local 
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Small Particles Large Particles 




1 2 3 4 1 2 3 

yld s y/d, 



FIG. 10: An image representation comparing the number den- 
sity distributions of 3D configurations for many different h. 
Black pixels represents a relative value of 1 and white repre- 
sent a relative value of 0. A gray scale is used to represent 
relative values between and 1. The pixel widths are 0.1d s 
horizontally and 0.2 vertically. 



packing fraction as a function of distance from the silo 
wall. Similar to our results, this local packing fraction 
showed fluctuations that decayed monotonically. In their 
paper they state a decay length of w 4d;; however, it 
appears that they drew this conclusion by estimating the 
value by eye. Applying Eqn. [2] to their data we find 
A on the order of di, close to the value found in our 
simulations. The experimental study by Seidler et al. 
reported on the local bond orientational order parameter 
which showed oscillation that decayed with distance from 
the wall. They reported a decay length of A ~ di using an 
exponential fit. The length scales from these two studies 
arc slightly larger than those found in our work. 



IV. MODEL 

Our results for <f>(h) can be understood with a simple 
model incorporating an effective boundary layer and a 
bulk like region. Figure [TT] shows a configuration of par- 
ticles confined between two plates and divided into three 
regions. Near each confining wall particles show strong 
layering and will assume a configuration much different 
than rep. For this model, these particles will be replaced 
by an effective "boundary layer" with packing fraction 
4>i- Layering will persist far into the bulk, but after some 
distance, indicated in the figure by 5L, the local pack- 
ing of particles will be near that of rep. Particles in the 
bulk region will be approximated to be in a rep config- 
uration with packing fraction <f> rcp . Using this simple 
model, 4> can be approximated by the weighted average 




FIG. 11: (Color online). Illustration of our model for <j>(h). 
Our model breaks a configuration with confining width h into 
three regions. The boundary layers are approximated to have 
a packing fraction <f)i and persist a distance 5L into the sam- 
ple, and the middle "bulk" region is approximated to have a 
packing fraction (j> rC p. 



different values of the parameters depending on the di- 
mension.) Reducing this equation further we obtain 



C 



(5) 



C = 0.233 



where C = 2SL(<j) rcp — <j>i). Note that this is the same 
form for <p{h) obtained from considering a 1st order cor- 
rection in terms of 1/h. 

Equation [5] only contains two fitting parameters, (f> rC p 
and C . 4> rC p is the packing fraction of an infinite uncon- 
fined configuration, and C approximates the difference in 
from 4> rcp for the first layers contained within a distance 
5L from the confining wall. The data in both Fig. [3] and 
Fig. |] are fitted to Eqn. [5| The fits are shown as the red 
lines (dark gray) in the earlier figures, and also in Fig.[l2| 
where the data are plotted as functions of 1/h to better 
illustrate the success of this model. The fits give for 2D 
0.844 and C = 0.317 and for 3D (f> rcp = 0.646 and 
Both fits give values for <j> rcp that are slightly 
larger, but not by much, than </> rcp reported earlier in 
the paper. In Fig. [12] is can be seen that the packing 
fraction for large 1/h dip significantly below the fitting 
line, due to the fluctuations in <fi(h) for small h; this is 
responsible for the over estimate in (f> rcp . When the data 
for both curves are fitted for h > 8 (1/h < 0.125) the 
actual values for 4> rcp are obtained. 

To provide further credence to our model we also per- 
form 2D rep simulations with a fixed circular boundary. 
Our simulations were carried out with different confining 
widths ranging from h 10 — 40, where h is the diame- 
ter of the circular boundary normalized by d s . Figure [13] 
shows a plot of 4>(h) for a circular fixed boundary. As 
before with a flat boundary condition, we see that <f> in- 
creases to an asymptotic limit. Adapting our model to a 
curved boundary with radius R = h/2, the weighted av- 
erage becomes <j>(R) = * [R ~^ SL) 1 4>rc P + 7r( ^. fl t L) 4>i- 
This expression can be simplified as <fi = 4> rcp — 2C/h, 
using C — 2SL(cj) r cp — <fti) as before, and dropping a 
term that is second order in SL/h. We show this curve 
as the red line (light gray) curve in Fig. 



13 using the 



h-2SL , 



-^j^4>i (in either 2D or 3D, with of course 



same values from our prior 2D fit with non-curved walls 
(^rcp = 0.844, C = 0.317), and find good agreement 
A direct fit to the circular data gives 



J rcp 

with the data 
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0.70 




0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 
1/h 

FIG. 12: (Color online). The upper black curve is a plot of 
<j>(l/h) for 2D configurations, and the red (dark gray) line 
going through the curve, is a fit from our model. Likewise, 
the lower black curve is a plot of <f>(l/h) for 3D configurations 
with the red (dark gray) line going through the curve being 
another fit from our model. 
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FIG. 13: (Color online). The black curve is the packing frac- 
tion dependence of random close packed 2D disks enclosed 
within a circular boundary. The red (light gray) curve is from 
our model using the parameters appropriate to Fig. 12 (2D 
case). The image at the lower right is an rep configuration 
confined in a circular boundary with h — 21. Small particles 
are rendered as green (medium gray) and large particles are 
rendered as blue (dark gray). 



gating simulated rep configurations confined between two 
walls in 2D and 3D. We find that confinement lowers the 
packing fraction, and induces heterogeneity in particle 
density where particles layer in bands near the wall. The 
structure of the local packing decays from a more or- 
dered packing near the wall to a less ordered packing in 
the bulk. All measures of local order and local density de- 
cay rapidly to their bulk values with characteristic length 
scales on the order of particle diameters. Thus, the influ- 
ence of the walls is rapidly forgotten in the interior of the 
sample, with confinement having the most notable effects 
when the confining dimension is quite small, perhaps less 
than 10 particle diameters across. 

These findings have implications for experiments inves- 
tigating the dynamics of densely packed confined systems 
(i.e. colloidal suspensions or granular materials) . For ex- 
ample, our work shows that for small h the packing frac- 
tion has significant variations at small h (mostly clearly 
seen in 2D, for example Fig. [3]). For dense particulate 

, flow is already difficult. By 



< 



J rcp 



suspensions with 
choosing a value of h with a local maximum in rcp (/i), 
a suspension may be better able to flow, as there will be 
more free volume available. Likewise, a poor choice of 
h may lead to poor packing and enhanced clogging. A 
microfluidic system with a tunable size h may be able to 
vary the flow properties significantly with small changes 
of h, but our work implies that control over h needs to 
be fairly careful to observe these effects. Of course, these 
effects will be obscured by polydispersity in many sys- 
tems of practical interest; however, our work certainly 
has implications for microfluidic flows of these sorts of 
materials, once the minimum length scales approach the 
mean particle size. 

Our work has additional implications for experiments 
on confined glasses [111 [SUl Ell E21 E31 EH EU [^1 [57] . As 
mentioned in the introduction, confinement changes the 
properties of glassy samples, but it is unclear if this is 
due to finite size effects or due to interfacial influences 
from the confining boundaries [58 . Our results show 
that dense packings have significant structural changes 
near the flat walls, suggesting that indeed interfacial in- 
fluences on materials can be quite strong at very short 
distances, assuming that the structural changes couple 
with dynamical behavior. 



4>rcp — 0.843 and C — 0.301; the slightly smaller value 
of C suggests that particles pack more efficiently near a 
curved boundary than a flat boundary. However, within 
the uncertainty of our limited data in the curved geome- 
try, it is not clear if the difference in C is significant. Acknowledgments 
V. CONCLUSION 



In this paper, we have shown how a confining boundary 
alters the structure of random close packing by investi- 



This work was supported by the National Science 
Foundation under Grant No. DMR-0804174. 



10 



[1] F. Zamponi, Nature 453, 606 (2008). 

[2] J. C. Hales, Ann. Math. 162, 1065 (2005). 

[3] G. D. Scott and D. M. Kilgour, J. Phys. D Appl. Phys. 

2, 863 (1969). 
[4] J. G. Berryman, Phys. Rev. A 27, 1053 (1983). 
[5] G. Y. Onoda and E. G. Liniger, Phys. Rev. Lett. 64, 

2727 (1990). 

[6] J. D. Bernal and J. Mason, Nature 188, 910 (1960). 
[7] C. S. O'Hern, S. A. Langer, A. J. Liu, and S. R. Nagel, 

Phys. Rev. Lett. 88, 075507 (2002). 
[8] N. Xu, J. Blawzdziewicz, and C. S. O'Hern, Phys Rev E 

71, 061306 (2005). 
[9] S. Torquato, T. M. Truskett, and P. G. Debenedetti, 

Phys. Rev. Lett. 84, 2064 (2000). 
[10] W. S. Jodrey and E. M. Tory, Phys. Rev. A 32, 2347 

(1985). 

[11] J. Tobochnik and P. M. Chapin, J. Chem. Phys. 88, 5824 
(1988). 

[12] A. R. Kansal, S. Torquato, and F. H. Stillinger, J. Chem. 

Phys. 117, 8212 (2002). 
[13] R. Al-Raoush and M. Alsalch, Powder Technol. 176, 47 

(2007) . 

[14] K. Lochmann, L. Oger, and D. Stoyan, Solid State Sci. 
8, 1397 (2006). 

[15] P. Richard, L. Oger, J. P. Troadec, and A. Gervois, Eur. 
Phys. J. E 6, 295 (2001). 

[16] H. J. H. Brouwers, Phys. Rev. E 74, 031309 (2006). 

[17] A. Doncv, I. Cisse, D. Sachs, E. A. Variano, F. H. Still- 
inger, R. Connelly, S. Torquato, and P. M. Chaikin, Sci- 
ence 303, 990 (2004). 

[18] T. Okubo and T. Odagaki, J Phys.-Condens. Mat. 16, 
6651 (2004). 

[19] K. Desmond and S. V. Franklin, Phys. Rev. E 73, 031306 
(2006). 

[20] Y. Jiao, F. H. Stillinger, and S. Torquato, Phys. Rev. 

Lett. 100, 245504 (2008). 
[21] S. F. Edwards, Granular Matter (Springer- Verlag, 1994). 
[22] C. Radin, J. Stat. Phys. 131, 567 (2008). 
[23] P. Jalali and M. Li, J. Chem. Phys. 120, 1138 (2004). 
[24] C. Song, P. Wang, and H. A. Makse, Nature 453, 629 

(2008) . 

[25] R. Pal, Polym. Eng. Sci. 48, 1250 (2008). 

[26] G. Lois, J. Blawzdziewicz, and C. S. O'Hern, 

arXiv:0809.1044 (2008). 
[27] R. Zallen, Physics of Amorphous Solids (Wiley, 1983). 
[28] M. Alcoutlabi and G. B. McKenna, J Phys.-Condens. 

Mat. 17, R461 (2005). 
[29] P. Gallo, M. Rovere, and E. Spohr, J. Chem. Phys. 113, 

11324 (2000). 
[30] S. Granick, Science 253, 1374 (1991). 
[31] J. R. Henderson, Mol. Phys. pp. 2345-2352 (2007). 
[32] J. Mittal, V. K. Shen, J. R. Errington, and T. M. 

Truskett, J. Chem. Phys. 127, 154513 (2007). 
[33] J. Mittal, J. R. Errington, and T. M. Truskett, J. Chem. 

Phys. 126, 244708 (2007). 
[34] J. Mittal, T. M. Truskett, J. R. Errington, and G. Hum- 
mer, Phys. Rev. Lett. 100, 145901 (2008). 
[35] I. Cohen, T. G. Mason, and D. A. Weitz, Phys. Rev. Lett. 

93, 046001 (2004). 
[36] K. To, P.-Y. Lai, and H. K. Pak, Physica A 315, 174 

(2002). 



[37] I. Zuriguel, L. A. Pugnaloni, A. Garcimartm, and 

D. Maza, Phys. Rev. E 68, 030301 (2003). 

[38] S. Redner and S. Datta, Phys. Rev. Lett. 84, 6018 (2000). 
[39] S. B. Fuller, E. J. Wilhelm, and J. M. Jacobson, J. Mi- 

croelectromech S. 11, 54 (2002). 
[40] K. Sharp and R. Adrian, Microfluid Nanofluid 1, 376 

(2005). 

[41] G. M. Whitesides, Nature 442, 368 (2006). 

[42] L. Vanel, P. Claudin, Bouchaud, M. E. Cates, 

E. Clement, and J. P. Wittmer, Phys. Rev. Lett. 84, 
1439 (2000). 

[43] U. Marconi, Physica A 280, 279 (2000). 
[44] P. Claudin and J.-P. Bouchaud, Phys. Rev. Lett. 78, 231 
(1997). 

[45] J. W. Landry, G. S. Grest, and J. Plimpton, Powder 
Technol. 139, 223 (2004). 

[46] J. W. Landry, G. S. Grest, L. E. Silbert, and S. J. Plimp- 
ton, Phys. Rev. E 67, 041303 (2003). 

[47] G. T. Seidler, G. Martinez, L. H. Seeley, K. H. Kim, 
E. A. Behne, S. Zaranek, B. D. Chapman, S. M. Heald, 
and D. L. Brewe, Phys. Rev. E 62, 8175 (2000). 

[48] L.-W. Teng, P.-S. Tu, and I. Lin, Phys. Rev. Lett. 90, 
245004 (2003). 

[49] C. Murray, MRS Bulletin 23, 33 (1998). 

[50] Z. T. Nemeth and H. Lowen, Phys. Rev. E 59, 6824 
(1999). 

[51] C. R. Nugent, K. V. Edmond, H. N. Patel, and E. R. 

Weeks, Phys. Rev. Lett. 99, 025702 (2007). 
[52] D. Morineau, Y. Xia, and C. A. Simionesco, J. Chem. 

Phys. 117, 8966 (2002). 
[53] J. Schiiller, Y. Mel'nichenko, R. Richert, and E. W. Fis- 
cher, Phys. Rev. Lett. 73, 2224 (1994). 
[54] C. L. Jackson and G. B. McKenna, J. Non-Cryst. Solids. 

131-133, 221 (1991). 
[55] K. Kim and R. Yamamoto, Phys. Rev. E 61, R41 (2000). 
[56] P. A. Thompson, G. S. Grest, and M. O. Robbins, Phys. 

Rev. Lett. 68, 3448 (1992). 
[57] P. Scheidler, W. Kob, K. Binder, and G. Parisi, Philos. 

Mag. A 82, 283 (2002). 
[58] F. He, L. M. Wang, and R. Richert, Eur. Phys. J. Special 

Topics 141, 3 (2007). 
[59] A. S. Clarke and J. D. Wiley, Phys. Rev. B 35, 7350 

(1987). 

[60] J. Nocedal and S. J. Wright, Numerical Optimization 

(Springer, 1999). 
[61] A. Donev, S. Torquato, F. H. Stillinger, and R. Connelly, 

Phys. Rev. E 70, 043301 (2004). 
[62] C. S. O'Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, 

Phys. Rev. E 70, 043302 (2004). 
[63] A. H. Marcus and S. A. Rice, Phys. Rev. Lett. 77, 2577 

(1996). 

[64] P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. 

Rev. B 28, 784 (1983). 
[65] U. Gasser, E. R. Weeks, A. Schofield, P. N. Pusey, and 

D. A. Weitz, Science 292, 258 (2001). 



